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Abstract. Arterial Spin Labelling (ASL) functional Magnetic Reso¬ 
nance Imaging (fMRI) data provides a quantitative measure of blood 
perfusion, that can be correlated to neuronal activation. In contrast to 
BOLD measure, it is a direct measure of cerebral blood flow. However, 
ASL data has a lower SNR and resolution so that the recovery of the per¬ 
fusion response of interest suffers from the contamination by a stronger 
hemodynamic component in the ASL signal. In this work we consider a 
model of both hemodynamic and perfusion components within the ASL 
signal. A physiological link between these two components is analyzed 
and used for a more accurate estimation of the perfusion response func¬ 
tion in particular in the usual ASL low SNR conditions. 


1 Introduction 

Arterial Spin Labelling (ASL) [1] provides a direct measure of cerebral blood 
flow (CBF), overcoming one of the most important limitations of Blood Oxygen 
Level Dependent (BOLD) signal [2]: BOLD contrast cannot quantify cerebral 
perfusion. In contrast to BOLD, ASL is able to provide a measure of baseline 
CBF as well as quantitative CBF signal changes in response to stimuli pre¬ 
sented to any volunteer in the scanner during an experimental paradigm. Hence, 
ASL enables the comparison of CBF changes between experiments and subjects 
(healthy vs patients) making its application to clinics feasible. In addition, ASL 
signal localization is closer to neural activity. ASL has already been used in clin¬ 
ics in steady-state for instance for probing CBF discrepancy in pathologies like 
Alzheimer’s disease and stroke, but its use in the functional MRI context has 
been limited so far. Despite ASL advantages, its main limitation lies in its low 
Signal-to-Noise Ratio (SNR), which, together with its low temporal and spatial 
resolutions, makes the analysis of such data more challenging. 

According to m, ASL signal has been typically analyzed with a general 
linear model (GLM) approach, accounting for a BOLD component mixed with 
the perfusion component. In such a setting both the hemodynamic response 
function (HRF or BRF for BOLD response function) and perfusion response 
function (PRF) are assumed to be the same and to fit the canonical BRF shape. 
In contrast, an adaptation of the Joint-Detection estimation (JDE) framework [5] 
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to ASL data has been proposed in m to separately estimate BRF and PRF 
shapes, and implicitly consider the control/label effect which, as stated in [4], 
increases the sensitivity of the analysis compared to differencing approaches. Al¬ 
though this JDE extension provides a good estimate of the BRF, the PRF esti¬ 
mation remains much more difficult because of the noisier nature of the perfusion 
component within the ASL signal. In the past decade, physiological models have 
been described to explain the physiological changes caused by neural activity. 
In [HU], neural coupling, which maps neural activity to ensuing CBF, and the 
Balloon models which relates CBF to BOLD signal, have been introduced. These 
models describe the process from neural activation to the BOLD measure, and 
the impact of neural activation on other physiological parameters. 

Here, we propose to rely on these physiological models to derive a tractable 
linear link between perfusion and BOLD components within the ASL signal and 
to exploit this link as a prior knowledge for the accurate and reliable recovery of 
the PRF shape in functional ASL data analysis. This way, we refine the separate 
estimation of the response functions in m by taking physiological information 
into consideration. The structure of this paper goes as follows: the physiological 
model and its linearization to find the PRF/BRF link are presented in section]^ 
Starting then from the ASL JDE model described in section we extend the 
estimation framework to account for the physiological link in section]^ Einally, 
results on artificial and real data are presented and discussed in sections [5][^ 

2 A physiologically informed ASL/BOLD link 

Our goal is to derive an approximate physiologically informed relationship be¬ 
tween the perfusion and hemodynamic response functions so as to improve their 
estimation in a JDE framework m- We show in this section that, although 
this relationship is an imperfect link resulting from a linearization, it provides 
a good approximation and allows to capture important features such as a shift 
in time-to-peak from one response to another. Eor a physiologically validated 
model, we use the extended balloon model described below. 


2.1 The extended Balloon model 

The Balloon model was first proposed in uni to link neuronal and vascular pro¬ 
cesses by considering the capillary as a balloon that dilates under the effect of 
blood fiow variations. More specifically, the model describes how, after some 
stimulation, the local blood fiow fin{t) increases and leads to the subsequent 
augmentation of the local capillary volume This incoming blood is strongly 
oxygenated but only part of the oxygen is consumed. It follows a local decrease of 
the deoxyhemoglobin concentration ^(t) and therefore a BOLD signal variation. 
The Balloon model was then extended in [8| to include the effect of the neuronal 
activity u{t) on the variation of some auto-regulated fiow inducing signal '0(t) 
so as to eventually link neuronal to hemodynamic activity. The global physio¬ 
logical model corresponds then to a non-linear system with four state variables 
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Fig. 1. Effect of the physiological parameters on the BRF (left) and PRF (right) shapes. 
The parameters values proposed in [8] are used except for one parameter whose identity 
and value is modified as indicated in the plot. 


/in, ^ 1 ^} corresponding to normalized flow inducing signal, local blood flow, 
local capillary volume, and deoxyhemoglobin concentration. Their interactions 
over time are described by the following system of differential equations: 


dfinit) 
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with initial conditions '0(0) = 0, /in(0) = 1^(0) = ^(0) = 1. Lower case notation 
is used for normalized functions by convention. The system depends on 5 hemo¬ 
dynamic parameters: r^, r/ and are time constants respectively for signal 
decay/elimination, auto-regulatory feedback from blood flow and mean transit 
time, w reflects the ability of the vein to eject blood, and Eq is the oxygen ex¬ 
traction fraction. Another parameter 77 is the neuronal efficacy weighting term 
that models neuronal efficacy variability. 

Once the solution of the previous system is found, Buxton et al m proposed 
the following expression that links the BOLD response h{t) to the physiological 
quantities considering intra-vascular and extra-vascular components: 


h{t) = vo[h{i - m) + k 2 {i - ^) + ^ 3(1 - 


( 2 ) 


where ki , k 2 and ks are scanner-dependent constants and Vq is the resting blood 
volume fraction. According to [ 10 ], /ci = TEq, k 2 = 2 and ks = 2 Eq — 0.2 at a field 
strength of 1.5T and echo time TE = 40ms. 

The physiological parameters used are the ones proposed by Friston et al in 
[ 8 ]: Vb = 0.02, = 1.25, r/ = 2.5, = 1, re = 0.2, Eq = 0.8 and r] = 0.5. 

The BRF and PRF generated using these parameters with the physiological 
model are shown in Fig. [funder the label “Friston 00” (dashed line). The rest 
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of the curves show the effect of changing the physiological parameters: is a 
scaling factor and causes non-linearities above a certain value; controls the 
signal decay, which is more or less smooth; the auto-regulatory feedback rj 
regulates the undershoot; the transit time expands or contracts the signal 
in time; the windkessel parameter w models the initial dip and the response 
magnitude; the oxygen extraction Eq impacts the response scale. After analysing 
the behaviour of the model when varying the parameters values, the impact of 
each parameter was investigated and we concluded that the values proposed in 
[8] seemed reasonable. 


2.2 Physiological linear relationship between response functions 

From the system of equations previously defined, we derive an approximate 
relationship between the PRF, namely g{t)^ and the BRF, which is given by 
h{t) when u{t) is an impulse function. Both BRF and PRF are percent signal 
changes, and we consider g{t) = fin{t) — 1, as fin{t) is the normalized perfusion, 
with initial value 1. Therefore the state variables are {'0,^', 1 — iz, 1 — 

In the following we will drop the time index t and consider functions h, -0, etc. 
in their discretized vector form. We can obtain a simple relationship between h 
and g by linearizing the system of equations. Equation <§ can first be linearized 
into: 


h = Vo[{h + k2){l - 0 + (fe - ^2)(1 - . (3) 

We then linearize the system 0 around the resting point {^p, g,l — 1^,1 — = 

{ 0 , 0 , 0 , 0 } as in m- From this linearization, denoting by V the first order dif¬ 
ferential operator and I the identity matrix, we get: 

'V{g} = -Ip 

, (P+^){1-*.} 

where 7 = ^ ^ follows a linear link between h and g 

that we write diS g = f2h where: 

= Vq-i (-{ki + k2hB + (fci + (5) 

V WtI Tm ) 

with A = (v -\- —— ) and B = (v ^ -) ( 6 ) 

V J \ TmJ 

Using values of physiological constants as proposed in El, Fig. I shows the 
BRF and PRF results that we get {hun^ gun) by applying the linear operator to 
physiologically generated PRF {gphysio) or BRF {hphysio)- hii„ = S2~^gphysio or 
giin = ^hphysio compared to these physiologically generated hphysio and gphysio 
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Fig. 2. Physiological responses generated with the physiological model, using param¬ 
eters proposed in [8]: neural activity ip, physiological (hphysio or BRFp/i^/sio) and lin¬ 
earized (hiin or BRFiin) BRFs, physiological (gphysio or PRFp/i^/sio) and linearized 
(giin or PRFan) PRFs. 


functions, computed by using the physiological model differential equations. Note 
that, although time-to-peak (TTP) values are not exact, the linear operator 
maintains the shape of the functions and satisfyingly captures the main features 
of the two responses. We considered a finer temporal resolution than TR for f2 
and, besides this, there is no direct dependence on the TR. 

The derivation of this linear operator gives us a new tool for analyzing the 
ASL signal, although this link is subject to caution as linearity assumption is 
strong and this linearization induces approximation error. 

3 Bayesian hierarchical model for ASL data analysis 

The ASL JDE model described in m assumes a partitioned brain into several 
functional homogeneous parcels each of which gathers signals which share the 
same response shapes. In a given parcel 7^, the generative model for ASL time 
series, measured at times {tn)n=i:N where tn = nTR, N is the number of scans 
and TR the time of repetition, with M experimental conditions, reads V j e V, 
\V\ = J: 


M 

yj= 2 aJ-X'^h + c'J^WX'^g+ P£j + ajw + bj (7) 

(a) (6) (c) (d) (e) 

The signal is decomposed into (a) task-related BOLD and (b) perfusion compo¬ 
nents given by the first two terms respectively; (c) a drift component Pij already 
considered in the BOLD JDE [5]; (d) a perfusion baseline term ajW which com¬ 
pletes the modelling of the perfusion component; and (e) a noise term. 

ASL fMRI data consists in the consecutive and alternated acquisitions of 
control and magnetically tagged images. The tagged image embodies a perfusion 
component besides the BOLD one, which is present in the control image too. 
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The BOLD component is noisier compared to standard BOLD fMRI acquisition. 
The control/tag effect is implicit in the ASL JDE model with the use of matrix 
W. More specifically, we further describe each signal part below. 

(a) The BOLD component: h g represents the unknown BRF shape, 

with size F + 1 and constant within V. The magnitude of activation or BOLD 
response levels are a = {a^^j g F,m = 1 : M}. 

(b) The perfusion component: It represents the variation of the perfusion 

from the baseline when there is task-related activity, g g represents the 

unknown PRF shape, with size F + 1 and constant within V. The magnitude of 
activation or perfusion response levels are c= {c^^j GF,m = 1 :M}. W mod¬ 
els the control/tag effect in the perfusion component, and it is further explained 
below. 

(a-b) Considering At < TR the sampling period of h and whose temporal 
resolution is assumed to be the same, X = n = 1 : A^, / = 0 : F} is a bi¬ 

nary matrix that encodes the lagged onset stimuli. In |6l7j . BRF and PRF shapes 
follow prior Gaussian distributions h ~ Af{0,VhXh) and g ~ M{O^VgXg)^ with 
covariance matrices and Xg encoding a constraint on the second order deriva¬ 
tive so as to account for temporal smoothness. The BOLD (BRLs) and perfusion 
(PRLs) response levels (resp. a and c) are assumed to follow different spatial 
Gaussian mixture models but governed by common binary hidden Markov ran¬ 
dom fields {qj^^j G V} encoding voxels’ activation = 1,0 for activated, resp. 
non-activated) states for each experimental condition m. This way, BRLs and 
PRLs are independent conditionally to q: p(a,c|gr). An Ising model on q in¬ 
troduces spatial correlation as in [SE]. For further interest please refer to [5]. 
Univariate Gamma/Gaussian mixtures were used instead in [ 12 ] at the expense 
of computational cost. The introduction of spatial modelling through hidden 
Markov random fields gave an improved sensitivity/specificity compromise. 

(c) The drift term: It allows to account for a potential drift and any other 

nuisance effect (e.g. slow motion parameters). Matrix P = [pi,... ,po] of size 
N X O comprises the values of an orthonormal basis (ie., Ptp = lo). Vector 
ij = 0 = 1 : Oy contains the corresponding unknown regression coefficients 

for voxel j. The prior reads ij ~ A/’(0, v^Io)- 

(b-d) The control/tag vector w (V-dimensional): It encodes the difference 
in magnetization signs between control and tagged ASL volumes. Wt^ = 1/2 if 
tn is even (control volume) and Wt^ = — 1/2 otherwise (tagged volume), and 
W = diag(tt;) is the diagonal matrix with w as diagonal entries. 

(d) The perfusion baseline: It is encoded by aj at voxel j. The prior reads 

aj - A/'(0, Va). 

(e) The noise term: It is assumed white Gaussian with unknown variance 
bj ~A/’( 0 , vijIn)- 

Hyper-parameters 0. Non-informative Jeffrey priors are adopted for {vb^V£^ 
Voc} and proper conjugate priors are considered for the mixture parameters of 
BRLs (Oa) and PRLs (Oc)- 
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4 A physiologically informed 2-steps inference procedure 

The BOLD component is known to have a higher SNR than the perfusion com¬ 
ponent in the ASL signal, and can be estimated with a higher confidence. The 
link g = f2h that we derived between both components can then be used to 
inform the PRF from the BRF. Using this link the other way around may not 
be satisfying as it may result in a contamination of by a noisier g. 

This effect has been noticed in the implementation of a physiologically in¬ 
formed Bayesian procedure, considering the generative model Q, and the fol¬ 
lowing priors for the PRF and BRF h ~ N'{0^Vh^h) and g\h ~ N'{f2h^VgUg), 
with Eh = Eg = (Z\t)^(T> 2 ^ 2 )~^- L >2 is the truncated second order finite dif¬ 
ference matrix of size (F — 1) x (F — 1) that introduces temporal smoothness, as 
in m, and Vh and Vg are scalars that we set manually. As seen in Fig. I^Mid- 
die], this approach does not yield satisfying results, not only for the perfusion 
component, but also for the BOLD one, compared to the model presented in 

m- 

We therefore propose to exploit the described physiological link in a two-step 
procedure, in which we first identify hemodynamics properties (/i, dj^), and then 
use the linear operator f2 and the previously estimated hemodynamic proper¬ 
ties to recover the perfusion component (^, cj^). This way, we avoid an arising 
contaminating effect of g on the estimation of /i, as in the one-step approach in 
Fig. [^Middle]. Each step is based on a Gibbs sampling procedure as in m- 


4.1 Hemodynamics estimation step A4i 

In a first step Ati, our goal is to extract the hemodynamic components and the 
drift term from the ASL data. In the JDE framework 0- it amounts to initially 
considering the perfusion component as a generalized perfusion term, including 
perfusion baseline and event-related perfusion response. The generative model 
0 for ASL time series can be equivalently written, by grouping the perfusion 
terms involving W = diag{w)^ as 

M / M \ 

Vi = E + pej+wi 2 cfx^g+ajl +bj (8) 

m=l \m=l / 

where we consider ajW = Wajl. Note that the hemodynamics components 
BRE h and the drift term ij can be estimated first, by segregating them from a 
general perfusion term and a noise term. However, the perfusion component is 
considered in the residuals, so as to properly estimate its different contributions 
in a second step Al 2 - 

Given the estimated ^ ^ we then compute residuals 

containing the remaining perfusion component: 


M 




" =yj- E 


A 




-Ml 


m=l 


(9) 
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4.2 Perfusion response estimation step AI 2 

From the residuals of the first step we estimate the perfusion component. 

The remaining signal is, according to 0 , Vi = 1 ; J, 

M 

yM, ^ ^ 2 dpwX^g + ajw + bj (10) 

m=l 

In this step, we introduce a prior on gf, to account for the already described 
physiological relationship g = f2h: 

g\h^^ Af{ f2h^\VgEg), with Sg = Ip ■ (11) 

The significance of the 2-step approach is to first preprocess the data to 
subtract the hemodynamic component within the ASL signal, as well as the 
drift effect, and to focus in a second step on the analysis of the smaller perfusion 
effect. In [4], differencing methods were used to subtract components with no 
interest in the perfusion analysis and directly analyse the perfusion effect in the 
time series. In contrast to these methods, we expect to disentangle perfusion from 
BOLD components by identifying all the components contained in the signal, 
and to recover them more accurately. 

5 Simulation results 

The generative model for ASL time series in section has been used to gen¬ 
erate artificial ASL data. A low SNR has been considered, with TR = 1 s, 
mean ISI = 5.03 s, duration 25 s. A/" = 325 scans and two experimental con¬ 
ditions (M = 2) represented with 20 x 20-voxel binary activation label maps 
corresponding to BRL and PRL maps shown in Fig. For both conditions: 
{aY'lqj = 1) ~ A/’(2.2,0.3) and {c^\qj = 1) ~ A/'(0.48, 0 . 1 ). Parameters were cho¬ 
sen to simulate a typical low SNR ASL scenario, in which the perfusion compo¬ 
nent is much lower than the hemodynamics component. A drift ij ~ A/'(0, IOJ 4 ) 
and noise variance = 7 were considered. BRF and PRF shapes were simulated 
with the physiological model, using the physiological parameters used in [ 8 ]. 

In a low SNR context, the PRF estimate retrieved by the former approach 
developed in |6l7j is not physilogically relevant as shown in Fig. |4[(c), Top]. In 
the case of a physiologically informed Bayesian approach, considering a single- 
step solution as in Fig. [^Middle], the perfusion component estimation is worse 
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(a) (b) (c) 



(d) 







Fig. 4. Results on artificial data. Top row: non-physiological version. Middle row: 
physiological 1-step version. Bottom row: physiological 2-steps version. (a,d): esti¬ 
mated BRL and PRL effect size maps respectively. The ground-truth maps for the BRL 
and PRL are depicted in Fig[^ (b,c): BRF and PRF estimates, respectively, with their 
ground truth. 


than for the approach described in m and the BRF estimation is also degraded 
owing to the influence of the noisier perfusion component during the sampling. 
In contrast, the 2-steps method proposed here delivers a PRF estimate very close 
to the simulated ground truth (see Fig.[^(c), Bottom] with a BRF which is well 
estimated too. 

In Fig.[^ the robustness of both approaches with respect to the noise variance 
is studied, in terms of BRF and PRF recovery. The relative root-mean-square- 
error (rRMSE) is computed for the PRF and BRF estimates, i.e. rRMSE^ = 
110 — 11/110^^’^^®^ II where 0 g {h,g}. We observed that maintaining a good 

performance in the BRE estimation, we achieved a much better recovery of the 
PRE for noise variances larger than = 1. Therefore, with the introduction 
of the physiological link between BRE and PRE, we have improved the PRE 
estimation. 

6 Real data results 

Real ASL data were recorded during an experiment designed to map auditory 
and visual brain functions, which consisted of = 291 scans lasting TR = 3 
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Fig. 5. Relative RMSE for the BRF and PRF and the two JDE versions, wrt noise 
variance Vb ranging from 0.5 to 30. 


with TE = IS ms, FoV 192 mm, each yielding a 3-D volume composed of 
64 X 64 X 22 voxels (resolution of 3 x 3 x 3.5 mm^). The tagging scheme used 
was PICORE Q2T, with TIi = 700 ms, TI 2 = 1700 ms. The paradigm was 
a fast event-related design (mean I SI = 5.1 s) comprising sixty auditory and 
visual stimuli. Two regions of interest in the right temporal lobe, for the auditory 
cortex, and left occipital lobe, for the visual cortex, were defined manually. 

Fig. I^b-c) depicts the response estimates superimposed to the canonical 
shape which is in accordance with the BRF estimates for both methods. Indeed, 
we consider here an auditory region where the canonical version has been fit¬ 
ted. Accordingly, the BRL maps (Fig. [^a)) also look alike for both methods. 
However, PRF estimates significantly differ and the effect of the physiologically- 
inspired regularization yields a more plausible PRF shape for the 2-steps ap¬ 
proach compared with the non-physiological JDE version. Results on PRL maps 
(Fig. [^d)) confirm the improved sensitivity of detection for the proposed ap¬ 
proach. In the same way, in the visual cortex. Fig. [^f-g) shows the BRF and 
PRF estimates, giving a more plausible PRF shape for the 2-steps approach, too. 
For the detection results (Fig. |^h)), the 2-steps approach seems also to provide 
a much better sensitivity of detection. 

7 Discussion and conclusion 

Starting from non-linear systems of differential equations induced by physio¬ 
logical models of the neuro-vascular coupling, we derived a tractable linear op¬ 
erator linking the perfusion and BOLD responses. This operator showed good 
approximation performance and demonstrated its ability to capture both real¬ 
istic perfusion and BOLD components. In addition, this derived linear operator 
was easily incorporated in a JDE framework at no additional cost and with a 
significant improvement in PRF estimation, especially in critical low SNR sit¬ 
uations. As shown on simulated data, the PRF estimation has been improved 
while maintaining accurate BRF estimation. Real data results seem to confirm 
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the better performance of the proposed physiological approach compared to its 
competing alternative. In terms of validation, future work will be devoted to 
intensive validation on whole brain analysis and multiple subjects. 
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Auditory cortex results 
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Visual cortex results 
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Fig. 6. Comparison of the two JDE versions on real data in the auditory and visual 
cortex, (top row in auditory and visual cortex results): non-physiological version, 
(bottom row in auditory and visual cortex results): physiological 2-steps version. 
(a,e) and (d,h): estimated BRL and PRL effect size maps, respectively. (b,f) and 
(c,g): BRF and PRF estimates, respectively. The canonical BRF is depicted as a black 
dashed line, while PRF and BRF estimated are depicted in solid red and blue lines, 
respectively. 


























